!>   purpose: this program gives instructions on how to use fingreen3d.
!>
!>  license:
!>
!>    this routine is part of fingreen3d.
!>
!>    fingreen3d is a free software package: you can redistribute it and/or modify it
!>    under the terms of the gnu lesser general public license as published by the
!>    free software foundation, either version 3 of the license, or (at your option)
!>    any later version.
!>
!>    you should have received a copy of the gnu general public license (either v3
!>    or later), along with fingreen3d. if not, see <http://www.gnu.org/licenses/>.
!>
!>  author:
!>
!>    yingyi liu on jul.20, 2018
!>
!>  reference:
!>
!>    yingyi liu et al. a reliable open-source package for performance evaluation of
!>    floating renewable energy systems in coastal and offshore regions. energy conversion
!>    and management, 2018
program tutorialprog

    use fin_green_3d_m
    use, intrinsic :: iso_fortran_env, only: rk => real64
    implicit none

    integer i, j, tag, nrh
    integer, parameter :: nk = 200
    real(rk) r, zf, zp, v, h, wvn(nk)
    real(rk) x, y, pi, g
    real(rk), allocatable :: rhp(:)
    complex(rk), allocatable :: excv(:, :), grn(:, :)

    data g, pi/9.807_rk, 3.141592653589793_rk/

    write (*, *) 'fingreen3d tutorial'
    write (*, *) '===================='
    write (*, '(a)', advance='no') '计算示例 1 或 2：'
    read (*, *) i

    select case (i)
    case (1)
        open (1, file='data/test01/inputdata.txt', status='old')
    case (2)
        open (1, file='data/test02/inputdata.txt', status='old')
    case default
        write (*, *) '输入错误！'
        stop
    end select

    open (2, file='grn(1).txt', status='unknown')
    open (3, file='grn(2).txt', status='unknown')
    open (4, file='grn(3).txt', status='unknown')

    write (*, *)
    write (*, '(2x,a)') '----------------------------------------------------------------------------------'
    write (*, '(2x,a)') '  welcome to the use of the software package fingreen3d.'
    write (*, '(2x,a)') '  this is a tutorial example for teaching you how to call it from a user program. '
    write (*, '(24x,30a,10x)') 'author: yingyi liu, kyushu university, jul.20, 2018'
    write (*, '(2x,a)') '----------------------------------------------------------------------------------'

    write (2, '(1x,a6,5(2x,a18))') 'r/h', 're(g)_newman', 'im(g)_newman', 're(g)_present', 'im(g)_present', 'relative error'
    write (3, '(1x,a6,5(2x,a18))') 'r/h', 're(gr)_newman', 'im(gr)_newman', 're(gr)_present', 'im(gr)_present', 'relative error'
    write (4, '(1x,a6,5(2x,a18))') 'r/h', 're(gz)_newman', 'im(gz)_newman', 're(gz)_present', 'im(gz)_present', 'relative error'

    write (*, *)

!     read the input file

    read (1, *)
    read (1, *) zf, zp, v, h
    read (1, *) tag, nrh
    print *, 'echo input parameters:'
    write (*, '(2x,a6,f10.6,2x,a6,f10.6,2x,a6,f10.6,2x,a6,f10.6)') 'zf=', zf, 'zp=', zp, 'v=', v, 'h=', h
    write (*, '(4x,a6,i2,4x,a6,i4,4x,a16,f10.6,a4)') 'tag=', tag, 'nrh=', nrh, 'wave frequency=', dsqrt(g*v)/(2._rk*pi), 'hz'

    allocate (rhp(nrh), excv(nrh, 3), grn(nrh, 3))
    call dispersion(wvn, nk, dsqrt(g*v), h)

    do i = 1, 12
        read (1, *)
    end do

    do j = 1, 3
        read (1, *)
        read (1, *)
        do i = 1, nrh
            read (1, *) rhp(i), x, y
            excv(i, j) = dcmplx(x, y)
        end do
    end do

!     calculate the green function by calling fingreen3d

    do i = 1, nrh

        r = rhp(i)*h

        call fingreen3d(r, zf, zp, v, wvn, nk, h, grn(i, :), tag)

        write (2, 1020) r/h, excv(i, 1), grn(i, 1), cdabs(excv(i, 1) - grn(i, 1))/cdabs(excv(i, 1))
        write (3, 1020) r/h, excv(i, 2), grn(i, 2), cdabs(excv(i, 2) - grn(i, 2))/cdabs(excv(i, 2))
        write (4, 1020) r/h, excv(i, 3), grn(i, 3), cdabs(excv(i, 3) - grn(i, 3))/cdabs(excv(i, 3))

    end do

    do j = 1, 3
        print *
        print *
        if (j == 1) then
            write (*, '(1x,a5,5(2x,a14))') 'r/h', 're(g)_newman', 'im(g)_newman', 're(g)_present', 'im(g)_present', 'relative error'
        else if (j == 2) then
        write (*, '(1x,a5,5(2x,a14))') 'r/h', 're(gr)_newman', 'im(gr)_newman', 're(gr)_present', 'im(gr)_present', 'relative error'
        else if (j == 3) then
        write (*, '(1x,a5,5(2x,a14))') 'r/h', 're(gz)_newman', 'im(gz)_newman', 're(gz)_present', 'im(gz)_present', 'relative error'
        end if
        do i = 1, nrh
            write (*, 1050) rhp(i), excv(i, j), grn(i, j), cdabs(excv(i, j) - grn(i, j))/cdabs(excv(i, j))
        end do
    end do

    deallocate (rhp, excv, grn)

    print *
    print *, 'the test program runs successfully.'

1020 format(f8.4, 5e20.12)
1050 format(2x, f8.4, 5e14.4)

contains

    !>   purpose: this subroutine computes roots of the water-wave dispersion equation
    !>             in finite water depth, by using a higher-order iterative method
    !>
    !>  license:
    !>
    !>    this routine is part of fingreen3d.
    !>
    !>    fingreen3d is a free software package: you can redistribute it and/or modify it
    !>    under the terms of the gnu lesser general public license as published by the
    !>    free software foundation, either version 3 of the license, or (at your option)
    !>    any later version.
    !>
    !>    you should have received a copy of the gnu general public license (either v3
    !>    or later), along with fingreen3d. if not, see <http://www.gnu.org/licenses/>.
    !>
    !>  author:
    !>
    !>    yingyi liu on mar.23, 2017
    !>
    !>  reference:
    !>
    !>    yingyi liu et al. a reliable open-source package for performance evaluation of
    !>    floating renewable energy systems in coastal and offshore regions. energy conversion
    !>    and management, 2018
    !>
    !>    j.n. newman
    !>    numerical solutions of the water-wave dispersion relation
    !>    applied ocean research 12 (1990) 14-18
    !>
    !>  parameters:
    !>      input:   nrt --- integer, the number of roots required
    !>                w   --- real, wave angular frequency
    !>                h   --- real, water depth (h>0)
    !>      output:  wvn --- real, an array storing roots of the dispersion equation
    subroutine dispersion(wvn, nrt, w, h)
!
!   evaluation of the roots of the following equations
!   by higher-order iterative method
!   first root stored in wvn is from eq. (i)
!   the rest roots are from eq. (ii)
!   i) w*w/g = k tanh ( kh )
!   ii) -w*w/g = um tan ( umh )
!

        integer, intent(in) :: nrt
        real(rk), intent(in) :: w, h
        real(rk), intent(out) :: wvn(1:nrt)
        integer i, m
        real(rk) t, x, u, y, dnm, g, pi
        real(rk) fun, dfun, d2fun, trial, exx

        data g, pi/9.807_rk, 3.141592653589793_rk/

!------------------------------------------------------------------
! i. calculation of wave number (root of eq. (i))
!------------------------------------------------------------------
!
!   initialize iteration by an accurate chebyshev approximation
!   if y=x, use the approximation directly insteady of iteration
!   to avoid the singularity in the denomenator of the transcendental
!   function; otherwise, do the iterative procedure.
!
        x = w*w*h/g
        if (x > 0._rk .and. x <= 2._rk) then
            y = dsqrt(x)*(0.9994_rk + 0.1701_rk*x + 0.0305*x*x)
        else
            t = x*dexp(-2._rk*x)
            y = x + 2._rk*t - 6._rk*t*t
        end if

        if (dabs(y - x) < 1.e-10) then
            wvn(1) = x/h
        else
            m = 0
            exx = 1._rk
            do while (exx > 1.0d-10)
                trial = y
                dnm = trial*trial - x*x
                fun = dlog((trial + x)/(trial - x))/2._rk - trial
                dfun = -x/dnm - 1._rk
                d2fun = 2._rk*x*trial/(dnm*dnm)
                y = trial - fun/dfun*(1._rk + (fun/dfun)*(d2fun/dfun)/2._rk)
                exx = dabs(y - trial)
                m = m + 1
            end do
            wvn(1) = y/h
        end if

!------------------------------------------------------------------
! ii. calcultion of roots of eq. (ii), which characterizes
!     the evanescene modes in eigenfunction
!------------------------------------------------------------------
!
!   initialize iteration by a suitable starting approximation
!
        u = 3._rk*x/(7._rk + 3._rk*x)
        t = 0.0159_rk + 0.1032_rk*u + 4.3152_rk*u*u - 2.8768_rk*u*u*u
!
!   perform iterative procedure to find exact solution of um (m=1,..nrt-1)
!   of the transcendental equation eq. (ii)
!
        do i = 2, nrt
            m = 0
            exx = 1._rk
            do while (exx > 1.0d-10)
                trial = t
                y = (i - 1)*pi - trial
                dnm = y*y + x*x
                fun = atan2(x, y) - trial
                dfun = x/dnm - 1._rk
                d2fun = 2._rk*x*trial/(dnm*dnm)
                t = trial - fun/dfun*(1._rk + (fun/dfun)*(d2fun/dfun)/2._rk)
                exx = dabs(t - trial)
                m = m + 1
            end do
            y = (i - 1)*pi - t
            wvn(i) = y/h
            t = t - pi*x/(x*x + pi*i*(pi*(i - 1) - t))
        end do

    end subroutine dispersion

end program tutorialprog
